###################################################################################################
######### Script producing the tables and figures in the appendix to ##############################
######### Bureaucratic Turnover Under New Governments: ############################################
######### The moderating effect of organizational buffering #######################################
###################################################################################################
####preamble####
#packages
library(dplyr)
library(survival)
library(survminer)
library(stargazer)
library(lmtest)
library(coxed)
library(ggthemes)
library(viridis)
library(ggpubr)
library(splitstackshape)
library(reshape2)
library(car)
library(sandwich)
library(lubridate)
library(stringr)
library(fixest)
library(marginaleffects)
library(tidyr)
library(modelsummary)

#clean
rm(list = ls())

#setting directory for output
wd_out <- "figures/appendix/"

#### functions ####
# creating figures
print_plot <- function(plot=obj, name, width=12, height=8){
  ggsave(paste0(name, ".pdf"), plot = plot, width = width, height = height, 
         units = "in", path = paste0("./", wd_out))
  ggsave(paste0(name, ".png"), plot = plot, width = width, height = height, 
         units = "in", path = paste0("./", wd_out,"pngs"))
}

# creating tex tables
output <- function(name, format = "tex") {
  paste0(wd_out, name, ".", format)
}

#Function to calculate marginal effects for the stepwise models with 1 interaction term between change of government and a dummy variable
model_marginal_effects <- function(list, NO_FE = 0) {
  for (i in 1:length(list)) {
    for(conf in c(0.9, 0.95)) {
      temp <- marginaleffects::avg_slopes(list[[i]], 
                                          variables = names(list[[i]]$coefficients)[1+NO_FE], 
                                          by = names(list[[i]]$coefficients)[2+NO_FE],
                                          conf_level = conf) %>% 
        mutate(!!names(.[3]) := case_when(names(.[3]) == "bureaucratic_buffer" & .[3] == 1 ~ "Level-2", 
                                          names(.[3]) == "bureaucratic_buffer" & .[3] == 0 ~ "Level-1",
                                          names(.[3]) == "state_secretaries_NSDID" & .[3] == 1 ~ "Level-1 & Political buffer",
                                          names(.[3]) == "state_secretaries_NSDID" & .[3] == 0 ~ "Level-1 & No buffer",
                                          TRUE ~ NA),
               c_lvl = conf)
      difference <- marginaleffects::avg_slopes(list[[i]], 
                                                variables = names(list[[i]]$coefficients)[1+NO_FE], 
                                                by = names(list[[i]]$coefficients)[2+NO_FE],
                                                hypothesis = "sequential",
                                                conf_level = conf) %>% 
        select(p.value, estimate, std.error, conf.low, conf.high) %>% 
        mutate(!!names(list[[1]]$coefficients)[2+NO_FE] := "Difference-in-Differences",
               term = names(list[[i]]$coefficients)[1+NO_FE],
               c_lvl = conf)
      
      temp <- bind_rows(temp, difference) %>% 
        mutate(model = names(list[i]))
      
      if (i == 1 & conf == 0.9) {
        all_models <- temp
      } else {
        all_models <- full_join(all_models, temp)
      }
    }
  }
  return(all_models)
  
}

#Function to calculate marginal effects for more complex models, specified hypotheses and conditional variable
model_marginal_effects_mechanisms <- function(model, conditional_variables_str, hyp) {
  for(conf in c(0.9, 0.95)) {
    temp <- marginaleffects::avg_slopes(model, 
                                        variables = names(model$coefficients)[1], 
                                        by = conditional_variables_str,
                                        conf_level = conf) %>% 
      mutate(!!names(.[3]) := case_when(names(.[3]) == "bureaucratic_buffer" & .[3] == 1 ~ "Level-2", 
                                        names(.[3]) == "bureaucratic_buffer" & .[3] == 0 ~ "Level-1",
                                        names(.[3]) == "state_secretaries_NSDID" & .[3] == 1 ~ "Level-1 & Political buffer",
                                        names(.[3]) == "state_secretaries_NSDID" & .[3] == 0 ~ "Level-1 & No buffer", 
                                        names(.[3]) == "election_year" & .[3] == 1 ~ "Election year",
                                        names(.[3]) == "election_year" & .[3] == 0 ~ "Within-term",
                                        names(.[3]) == "bureaucratic_buffer_i" & .[3] == "Level-1 (Permanent Secretary)" ~ "Level-1 (Permanent Secretary)",
                                        names(.[3]) == "bureaucratic_buffer_i" & .[3] == "Level-2 (Director General)" ~ "Level-2 (Director General)",
                                        names(.[3]) == "bureaucratic_buffer_i" & .[3] == "Level-1 (Director General)" ~ "Level-1 (Director General)",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-1 NO POL" ~ "Level-1 & No Political buffer",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-2 and POL" ~ "Level-2 & Political buffer",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-1 and POL" ~ "Level-1 & Political buffer",
                                        names(.[3]) == "organizational_buffer2" & .[3] == "Level-2 NO POL" ~ "Level-2 & No Political buffer",
                                        .default = NA),
             c_lvl = conf)
    if (names(temp)[4] != "estimate") { 
      temp <- temp %>% 
        mutate(!!names(.[4]) := case_when(names(.[4]) == "bureaucratic_buffer" & .[4] == 1 ~ "Level-2", 
                                          names(.[4]) == "bureaucratic_buffer" & .[4] == 0 ~ "Level-1",
                                          names(.[4]) == "state_secretaries_NSDID" & .[4] == 1 ~ "Level-1 & Political buffer",
                                          names(.[4]) == "state_secretaries_NSDID" & .[4] == 0 ~ "Level-1 & No buffer", 
                                          names(.[4]) == "election_year" & .[4] == 1 ~ "Election year",
                                          names(.[4]) == "election_year" & .[4] == 0 ~ "Within-term",
                                          names(.[4]) == "bureaucratic_buffer_i" & .[4] == "Level-1 (Permanent Secretary)" ~ "Level-1 (Permanent Secretary)",
                                          names(.[4]) == "bureaucratic_buffer_i" & .[4] == "Level-2 (Director General)" ~ "Level-2 (Director General)",
                                          names(.[4]) == "bureaucratic_buffer_i" & .[4] == "Level-1 (Director General)" ~ "Level-1 (Director General)",
                                          .default = NA))
    }
    
    difference <- marginaleffects::avg_slopes(model, 
                                              variables = names(model$coefficients)[1], 
                                              by = conditional_variables_str,
                                              hypothesis = hyp,
                                              conf_level = conf) %>% 
      mutate(!!names(temp)[3] := "Difference-in-Differences",
             c_lvl = conf)
    
    if ("election_year" %in% conditional_variables_str) { 
      difference <- difference %>% 
        mutate(term = purrr::reduce2(c("1\\)", "0\\)"), c("Election year)", "Within-term)"), .init = term, str_replace_all))
    } else if ("state_secretaries_NSDID"  %in% conditional_variables_str) {
      difference <- difference %>% 
        mutate(term = purrr::reduce2(c("1\\)", "0\\)"), c("Political buffer)", "No Political buffer)"), .init = term, str_replace_all))
    }
    
    temp <- bind_rows(temp, difference)
    
    if (conf == 0.9) {
      all_models <- temp
    } else {
      all_models <- full_join(all_models, temp)
    }
  }
  return(all_models)
}

#Marginal effects function for Cox and Logistic models
cox_glm_marginal_effects_function <-  function(model, moderating_variabel, mod) { 
  
  combcoef <- c(1, 1, 0)*model$coefficients["GOV_turnover"]+
    c(0, 1, 1)*model$coefficients[str_c("GOV_turnover:", moderating_variabel)]
  
  if (mod == "Logistic") {
    vcov <- vcov(model, cluster ~ PersonID)
  } else { 
    vcov <- model$var 
    rownames(vcov) <- names(coefficients(model))
    colnames(vcov) <- names(coefficients(model))
  }
  
  combcoef.se <- c(sqrt(vcov["GOV_turnover","GOV_turnover"]),
                   sqrt(vcov["GOV_turnover","GOV_turnover"]+
                          (1^2)*vcov[str_c("GOV_turnover:", moderating_variabel), str_c("GOV_turnover:", moderating_variabel)]+
                          2*vcov["GOV_turnover", str_c("GOV_turnover:", moderating_variabel)]),
                   sqrt(vcov[str_c("GOV_turnover:", moderating_variabel), str_c("GOV_turnover:", moderating_variabel)]))
  
  group <- if (moderating_variabel == "bureaucratic_buffer") {
    group <- c("Level-1", "Level-2", "Difference-in-Differences")
  } else {
    group <- c("Level-1 & No buffer", "Level-1 & Political buffer", "Difference-in-Differences")
  }
  
  plotdata <- data.frame(name=moderating_variabel,
                         value=group,
                         estimate=exp(combcoef), 
                         conf.low=c(exp(combcoef+1.645*combcoef.se), exp(combcoef+1.96*combcoef.se)),
                         conf.high=c(exp(combcoef-1.645*combcoef.se), exp(combcoef-1.96*combcoef.se)),
                         c_lvl=c(rep("0.9", 3), rep("0.95", 3)),
                         term="Change of government",
                         model=mod)
  
  return(plotdata)
  
}


#Function to plot the marginal effects and the differences
marginal_effect_plot <- function(dataframe, facet = TRUE, models = FALSE, ORHR = FALSE, base_size = 16) {
  if (is.null(dataframe$name)) {
    dataframe <- dataframe %>% 
      pivot_longer(!!names(.[3]), values_drop_na = T)
  }
  
  if (facet == TRUE & models == FALSE & ORHR == FALSE) {
    plot <- ggplot(dataframe, aes(estimate, str_replace_all(term, output_names), xmin = conf.low, xmax = conf.high, 
                                  col = factor(value), linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
                                  group = value))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 2) +
      geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
      facet_wrap(~str_replace_all(name, output_names))+
      labs(x = "Marginal effect", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(0.4, 1.2))+
      scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  } else if (facet == TRUE & models == TRUE & ORHR == FALSE) {
    plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high, 
                                  col = model, linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
                                  group = model))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 2) +
      geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
      facet_wrap(~str_replace_all(name, output_names))+
      labs(x = "Marginal effect", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(0.4, 1.2))+
      scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  } else if (facet == TRUE & models == TRUE & ORHR == TRUE) {
    plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high, 
                                  col = model, linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
                                  group = model))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 2) +
      geom_vline(xintercept = 1, color = "darkred", linetype = 2)+
      facet_wrap(~str_replace_all(name, output_names))+
      labs(x = "Marginal effect (odds/hazard ratio)", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(1.2, 2.2))+
      scale_x_continuous(breaks = seq(0,6, by=0.5))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  } else {
    plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high,
                                  group = factor(value),
                                  linewidth = factor(c_lvl, levels = c("0.95", "0.9"))))+
      geom_linerange(position = position_dodge(0.4)) +
      geom_point(position = position_dodge(0.4), size = 3) +
      geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
      labs(x = "Marginal effect", y = NULL)+
      theme_base(base_size = base_size)+
      scale_color_viridis(discrete = T, option = "D", name = NULL)+
      scale_linewidth_discrete(range = c(1.4, 2.4))+
      scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
      guides(color = guide_legend(reverse=TRUE),
             linewidth = "none")+
      theme(legend.position = "top",
            legend.box = "vertical", plot.background=element_blank())
    print(plot) 
  }
}

#### Lists of names ####
output_names <- c("GOV_turnover_lead" = "Change of government (t+1)",
                  "GOV_turnover_lagged" = "Change of government (t-1)",
                  "GOV_turnover" = "Change of government",
                  "minister_turnover" = "Change of minister",
                  "bureaucratic_buffer_i" = "Bureaucratic buffer & position type",
                  "bureaucratic_buffer" = "Bureaucratic buffer",
                  "organizational_buffer2" = "Organizational buffer",
                  "state_secretaries_NSDID" = "Political buffer",
                  "state_secretaries" = "Political buffer",
                  "years_since_last_turnover_2" = "Time since change of government",
                  "years_since_last_turnover_2_categories" = "Time since change of government categorial",
                  "match_start_end_PMparty_flipped" = "Party incongruence",
                  "election_year" = "Election year",
                  "age_year" = "Age",
                  "gender" = "Gender",
                  "temporary_employment" = "Temporary appointment",
                  "pre_unionsoppløsning" = "Before-1906", 
                  "foreign_affairs" = "Ministry of Foreign Affairs",
                  "government_PMparty_share_posts_end_of_year" = "PM-party share of cabinet",
                  "years_since_DEP_start" = "Time since ministry created",
                  "dep_end" = "Ministry terminated",
                  "pol_flexible" = "Exposed to change of government",
                  "NSD_ID" = "Ministry Fixed Effects",
                  "decade" = "Decade",
                  "decade2" = "Decade",
                  "duration" = "Senior bureaucrat tenure in years",
                  "WWII" = "WWII",
                  "age_group_year2" = "Age Group",
                  "multiple_L1" = "Multiple Level-1 bureuacrats",
                  "dep.rad_turnover" = "Senior bureaucrat turnover")

moderators_names <- str_replace_all(sapply(c("Baseline",
                                             "government_PMparty_share_posts_end_of_year",
                                             "years_since_DEP_start", "dep_end",
                                             "multiple_L1", "age_group_year2"), 
                                           function(x){rep(x, 6)}), output_names)

DV_names <- lapply(c("Senior bureaucrat turnover", 
                     "Including ministry switching",
                     "Including ministry switching \n and promotion/demotion"), 
                   function(x){rep(x, 6)}) %>% unlist()

#### modelsummary fit statistics ####
gm <- tribble(~raw, ~clean, ~fmt,
              "FE: time", "FE: Time", 0,
              "FE: NSD_ID", "FE: Ministry", 0,
              "FE: decade", "FE: Decade", 0,
              "nobs", "N", 0,
              "r.squared", "R2", 2,
              "aic", "AIC", 2)

###################################
####preparing data for analysis####
###################################
#Load data
df <- readRDS(file = "data/TCS_NOR.RDS") 

#recoding
df$match_start_end_PMparty_flipped <- if_else(df$match_start_end_PMparty == 1, 0, 1)
df$election_year <- if_else(is.na(df$election_year), 0, df$election_year)
df$military <- ifelse(is.na(df$military), 0, df$military)

#### creating subsets of the data ####
raw <- df %>% filter(year >= 1884, military == 0) #removing military DGs and observations before parliamentarism

#pensioners are removed | and WWII
df <- raw %>%
  mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>% 
  filter(age_year <= 70, WWII == 0) 

#pensioners are removed 
df_med_WWII <- raw %>%
  mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>% 
  filter(age_year <= 70) 

#spells must be more than 2 years and pensioners are removed 
df2 <- df %>% 
  filter(less_than_two_years == 0) %>%
  mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>% 
  filter(age_year <= 70) 

#only Level-1 bureaucrats
df_level1 <- df %>% filter(bureaucratic_buffer == 0)

#### dataframe with all used variables
sub_df <- df_med_WWII %>% 
  select(start, stop, dep.rad_turnover, GOV_turnover, years_since_last_turnover_2, match_start_end_PMparty_flipped,
         election_year, age_year, gender, temporary_employment, pol_flexible, pre_unionsoppløsning, WWII, decade2,
         dep_end, years_since_DEP_start, government_PMparty_share_posts_end_of_year,
         foreign_affairs, WWI, post_WWII, NSD_ID, bureaucratic_buffer, year, 
         state_secretaries, proximity_to_GOV_turnover,proximity_to_GOV_turnover_2_categories, years_since_last_turnover_2_categories,
         GOV_turnover_any_change, GOV_turnover_new_PM, GOV_turnover_new_Parties,
         GOV_turnover_during_term, GOV_turnover_election, government_percent_parliament_end_of_year, 
         government_single_party_end_of_year, government_minority_end_of_year, government_left_right_gov_end_of_year,
         gov_turnover_during_spell, pension_age_propensity_age_year, pension_age_propensity_age_at_start,
         match_start_end_party, match_start_end_left_right, match_start_end_coalition_party,
         minister_turnover, political_career,
         GOV_turnover_lagged, GOV_turnover_lead, id_spells, duration, bureaucratic_buffer_i, 
         y_since_PS_in_min, decade_since_PS_in_min,rounded_y_since_PS_in_min, organizational_buffer, 
         organizational_buffer2, after_political_appointees_NSDID, after_state_secretaries_NSDID, after_poladv_NSDID,
         next_position_sector, education_main, education_level, post_SLS, pension_age_propensity_age_year,
         age_group_year2, multiple_L1, ministry_portfolio)

####Figure 9 - Types of government in Norway 1884-2024 ####
governments <- readRDS("data/Governments_NOR.RDS") %>% 
  mutate(government_type = case_when(regjeringstype == "Flertall ettparti" ~ "Single-Party Majority",
                                     regjeringstype == "Flertall minste vinnende koalisjon" |
                                       regjeringstype == "Flertall overtallig koalisjon"  ~ "Multi-Party Majority",
                                     regjeringstype == "Mindretall ettparti" ~ "Single-Party Minority",
                                     regjeringstype == "Mindretall koalisjon" ~ "Multi-Party Minority"))

#time labels
table <- governments %>% 
  group_by(government_type) %>% 
  summarise(time_with_government_type = sum(lengde_nesten_alle_endringer, na.rm = F),
            time_with_government_type_percentage = time_with_government_type / 
              (as.numeric( ymd("2024-12-31") - ymd("1884-06-26") ) ))

time_labels <- paste(table$government_type, 
                             paste0( paste0("(", scales::percent( as.numeric( table$time_with_government_type_percentage) ) ), ")" )
                             , sep = " ")

#plot
obj <- ggplot(governments, 
              aes(xmin = tiltrådt,
                  xmax = avgang,
                  ymin = 0,
                  ymax = percent,
                  fill = government_type))+
  geom_rect()+
  geom_hline(yintercept =0.5, linetype = 2, col = "red")+
  scale_x_date(expand = c(0.025, 0), date_breaks = "5 years", date_labels = "%Y", limits = ymd(c("1883-01-01", "2025-01-31") ) )+
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.10), labels = scales::percent)+
  scale_fill_viridis(discrete = T, labels = time_labels)+
  theme_base(base_size = 16)+
  labs(fill = NULL, y = "Percentage of seats in parliament", x = NULL) +
  theme(legend.position = "top",
        legend.box = "vertical",
        plot.background=element_blank(), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ 
  guides(fill=guide_legend(nrow = 3,
                           byrow = T,
                           title.position = "top"
  ))

print_plot(plot = obj, name = "governmentsovertime")

##Time with different governments
#minority
18.5+36.4
#single-party
23.8+36.4

##number of change in prime minister and party
sum(governments$statsminister_parti != lag(governments$statsminister_parti), na.rm = T)
# number of change after an election
sum(governments$lag_årsak_til_avgang[governments$statsminister_parti != lag(governments$statsminister_parti)] != "Valg", na.rm = T)
# Type of governments formed
table(governments$minority_majority[governments$statsminister_parti != lag(governments$statsminister_parti)])
# Type of governments formed
table(governments$koalisjon_ettparti[governments$statsminister_parti != lag(governments$statsminister_parti)])


####Figure 10 - The Percentage of ministries employing at least 1 state secretary and/or 1 permanent secretary over time ####
plotdata <- df %>% 
  distinct(year, NSD_ID, PS_in_dep, state_secretaries_NSDID, GOV_turnover) %>% 
  reframe(share_dep_rad = sum(PS_in_dep == 1)/n(),
          share_statesec = sum(state_secretaries_NSDID == 1)/n(),
          GOV_turnover_year = if_else(GOV_turnover == 1, year, NA),
          .by = year) %>% 
  pivot_longer(cols = share_dep_rad:share_statesec)

variance_in_positions <- ggplot(plotdata, aes(year, value, col = name))+
  geom_segment(aes(x = GOV_turnover_year, y = 0, yend = 0.03), col="black")+
  geom_line(size=1.3)+
  ylab("Percentage of Ministries with Position") +
  xlab(NULL) +
  theme_base(base_size = 26)+
  scale_color_viridis(discrete = T, option = "D", name = NULL, labels = c("Permanent Secretary", "State Secretary"))+
  scale_y_continuous(labels = scales::percent)+
  scale_x_continuous(breaks = seq(1884,2025, by=20))+
  guides(color = guide_legend(reverse=TRUE))+
  theme(legend.position = "top",
        legend.box = "vertical", plot.background=element_blank())

print_plot(variance_in_positions, name = "variance_in_positions")


####Table 3-5 - Descriptive Statistics ####
descriptive_stats <- sub_df %>% 
  group_by(id_spells) %>% 
  dplyr::select(duration, dep.rad_turnover, GOV_turnover, bureaucratic_buffer, state_secretaries,
                election_year, age_year, gender, temporary_employment, pol_flexible, foreign_affairs,
                pre_unionsoppløsning, government_PMparty_share_posts_end_of_year,
                years_since_DEP_start, dep_end) %>% 
  mutate(duration = round(as.numeric(duration/365.25))) %>%
  ungroup() %>% 
  dplyr::select(-id_spells)
names(descriptive_stats) <- str_replace_all(names(descriptive_stats), output_names)

#All
stargazer(as.data.frame(descriptive_stats),
          omit.summary.stat = c("p25", "p75"),
          title = "Descriptiv statistics: All senior bureaucrats",
          align=F,
          no.space=T, 
          model.numbers=T,
          label = c("Tab:desc1"),
          type = "latex",
          out = output(name = "descriptivstats_complete", format = "tex") 
)

##Level-1
highdata <- as.data.frame(descriptive_stats %>%
                            filter(`Bureaucratic buffer` == 1) %>% 
                            select(-`Bureaucratic buffer`))
stargazer(highdata,
          omit.summary.stat = c("p25", "p75"),
          title = "Descriptiv statistics: Level-1 Bureaucrats",
          align=F,
          no.space=T, 
          model.numbers=T,
          label = c("Tab:desc2"),
          type = "latex",
          out = output(name = "descriptivstats_highest", format = "tex")
) 

##Level-2
lowdata <- as.data.frame(descriptive_stats %>%
                           filter(`Bureaucratic buffer` == 0) %>% 
                           select(-`Bureaucratic buffer`) )
stargazer(lowdata,
          omit.summary.stat = c("p25", "p75"),
          title = "Descriptiv statistics: Level-2 Bureaucrats",
          align=F,
          no.space=T, 
          model.numbers=T,
          label = c("Tab:desc3"),
          type = "latex",
          out = output(name = "descriptivstats_lowest", format = "tex")
) 


####Figure 11 - The number of senior bureaucrats that exited their position at different intervals of tenure in years####
df$censoring <- if_else(df$Sensurert == "Høyre", "Right-censoring", "Exit")
df_med_WWII$censoring <- if_else(df_med_WWII$Sensurert == "Høyre", "Right-censoring", "Exit")

length <- ggplot(df_med_WWII %>% 
                   distinct(id_spells, duration, censoring), 
                 aes(duration/365.25, fill = censoring))+
  geom_bar(width = 0.9)+
  scale_x_continuous(limits = c(0, 41), breaks = c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40))+
  scale_y_continuous(limits = c(0,120))+
  theme_base(base_size = 16)+
  scale_fill_viridis(discrete = T, name = "")+
  labs(x = "Senior Bureaucrat Tenure in Years", y = "Number of Senior Bureaucrats") +
  theme(legend.position = "top",
        legend.box = "vertical", plot.background=element_blank())

print_plot(length, name = "length")

##length of spells at different points in time
df %>% 
  distinct(id_spells, duration, Sensurert, decade_start) %>% 
  filter(Sensurert == "ikke sensurert") %>% 
  mutate(periods = case_when(decade_start < 1940 ~ "Before 1940",
                             decade_start >= 1950 & decade_start < 1980 ~ "1950-1980",
                             decade_start >= 1980 & decade_start < 2010 ~ "1980-2010",
                             TRUE ~ "After 2010")) %>% 
  summarise(median = median(duration)/365.25,
            .by = c(periods))


####Figure 12 - Kaplan-Meier survival curves for senior bureaucrats in years with vs. without change of government, subset by bureaucratic and political buffering####
#Bureaucratic Buffering
km1 <- survfit(Surv(start, stop, dep.rad_turnover) ~
                 (GOV_turnover) + (bureaucratic_buffer), df, robust = T, id = id_spells)

plotdata1 <- ggsurvplot(km1)$data.survplot
plotdata1$start <- plotdata1$time-1
plotdata1 <- plotdata1 %>% 
  tidyr::pivot_longer(cols = c(time, start), values_to = "time", names_to = NULL) %>% 
  mutate(x_variabel = if_else(bureaucratic_buffer == 1, "Level-2", "Level-1"),
         grp_variabel = "Bureaucratic buffer") %>% 
  select(-bureaucratic_buffer)

#Political Buffering
km2 <- survfit(Surv(start, stop, dep.rad_turnover) ~
                   (GOV_turnover) + (state_secretaries_NSDID), df %>% filter(bureaucratic_buffer == 0), robust = T, id = id_spells)
plotdata2 <- ggsurvplot(km2)$data.survplot
plotdata2$start <- plotdata2$time-1
plotdata2 <- plotdata2 %>% 
  tidyr::pivot_longer(cols = c(time, start), values_to = "time", names_to = NULL) %>% 
  mutate(x_variabel = if_else(state_secretaries_NSDID == 1, "Level-1 & Political buffer", "Level-1 & No buffer"),
         grp_variabel = "Political buffer") %>% 
  select(-state_secretaries_NSDID)

#Combining data
plotdata <- full_join(plotdata1, plotdata2)

#Ploting
obj <- ggplot(plotdata, aes(time, surv, 
                            ymax = upper, ymin = lower,
                            fill = GOV_turnover, col = GOV_turnover))+
  geom_line(linetype = 1, linewidth = 1)+
  ggh4x::facet_nested_wrap(vars(grp_variabel, x_variabel), nrow = 2, ncol = 2)+
  scale_x_continuous(limits = c(0, 41), breaks = seq(0, 40, by = 5))+
  scale_y_continuous(limits = c(0,1), labels = scales::percent)+
  theme_base(base_size = 16)+
  scale_color_viridis(discrete = T, name = "", labels = c("Within-term", "Change of government"), direction = -1)+
  scale_fill_viridis(discrete = T, name = "", labels = c("Within-term", "Change of government"), direction = -1)+
  labs(x = "Senior Bureaucrat Tenure in Years", y = "") +
  theme(legend.position = "top",
        legend.box = "vertical", plot.background=element_blank())

print_plot(obj, name = "kapmeierplots")


#### Figure 13 - Hazard/Odds ratios with 95% (small stroke width) and 90% (big stroke width) confidence bands for the marginal effect of change of government when estimating the effect with alternative estimators.  ####
## Bureaucratic buffer 
#Cox
cox_no_FE <- coxph(Surv(start, stop, dep.rad_turnover) ~ GOV_turnover * bureaucratic_buffer + 
                     election_year + age_year + gender + temporary_employment + pre_unionsoppløsning + as.factor(decade) + as.factor(NSD_ID),
                   cluster = PersonID, 
                   df)
modelsummary(cox_no_FE, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("cox_no_FE", "html"), 
             coef_rename = output_names,
             coef_omit = "decade|NSD_ID",
             exponentiate = TRUE, 
             gof_map = gm,
             add_rows = data.frame(rbind(cbind("FE: Time", " "),
                                         cbind("FE: Ministry", "X"),
                                         cbind("FE: Decade", "X"),
                                         cbind("Number of senior bureaucrats", length(unique(df$id_spells[is.finite(df$bureaucratic_buffer)]))),
                                         cbind("Number of events", cox_no_FE$nevent))),
             title = "Cox regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Hazard Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

#Logistic
GLM_no_FE <- feglm(dep.rad_turnover ~ GOV_turnover * bureaucratic_buffer + time + time_squared + time_cubed +
                     csw0(election_year, age_year, gender, temporary_employment, pre_unionsoppløsning) | decade + NSD_ID,
                   cluster = "PersonID", 
                   family = "binomial",
                   df)
modelsummary(GLM_no_FE, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("GLM_no_FE", "html"), 
             coef_rename = output_names,
             exponentiate = TRUE,
             col.names = NULL,
             gof_map = gm,
             title = "Step-wise Logistic regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Odds Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

## Political buffer
cox_no_FE_polbuff <- coxph(Surv(start, stop, dep.rad_turnover) ~ GOV_turnover * state_secretaries_NSDID + 
                             election_year + age_year + gender + temporary_employment + pre_unionsoppløsning + as.factor(decade) + as.factor(NSD_ID),
                           cluster = PersonID, 
                           df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
modelsummary(cox_no_FE_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("cox_no_FE_polbuff", "html"), 
             coef_rename = output_names,
             coef_omit = "decade|NSD_ID",
             exponentiate = TRUE,
             gof_map = gm,
             add_rows = data.frame(rbind(cbind("FE: Time", " "),
                                         cbind("FE: Ministry", "X"),
                                         cbind("FE: Decade", "X"),
                                         cbind("Number of senior bureaucrats", length(unique(df$id_spells[is.finite(df$bureaucratic_buffer) & df$bureaucratic_buffer_i != "Level-2 (Director General)"]))),
                                         cbind("Number of events", cox_no_FE_polbuff$nevent))),
             title = "Cox regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Hazard Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

GLM_no_FE_polbuff <- feglm(dep.rad_turnover ~ GOV_turnover * state_secretaries_NSDID + time + time_squared + time_cubed +
                             csw0(election_year, age_year, gender, temporary_employment, pre_unionsoppløsning) | decade + NSD_ID,
                           cluster = "PersonID", 
                           family = "binomial",
                           df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
modelsummary(GLM_no_FE_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("GLM_no_FE_polbuff", "html"), 
             coef_rename = output_names,
             exponentiate = TRUE,
             col.names = NULL,
             gof_map = gm,
             title = "Step-wise Logistic regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Odds Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

plotdata <- bind_rows(
  cox_glm_marginal_effects_function(cox_no_FE, "bureaucratic_buffer", mod = "Cox"),
  cox_glm_marginal_effects_function(cox_no_FE_polbuff, "state_secretaries_NSDID", mod = "Cox"),
  cox_glm_marginal_effects_function(GLM_no_FE[[6]], "bureaucratic_buffer", mod = "Logistic"),
  cox_glm_marginal_effects_function(GLM_no_FE_polbuff[[6]], "state_secretaries_NSDID", mod = "Logistic"))

print_plot(marginal_effect_plot(plotdata, models = T, ORHR = T), name = "Coefplot_logistic_cox")


#### Figure 14 - Marginal effects with 95\% (small stroke width) and 90\% (big stroke width) confidence bands for the marginal effect of change of government on alternative specifications of the dependent variable: senior bureaucrat turnover ####
## Bureaucratic buffer 
OLS_alternative_DVs <- feols(c(dep.rad_turnover, dep.rad_turnover2, dep.rad_turnover3) ~ GOV_turnover * bureaucratic_buffer +
                               election_year + age_year + gender + temporary_employment + pre_unionsoppløsning | time+NSD_ID+decade, 
                             cluster = "PersonID", df)
names(OLS_alternative_DVs) <- c("Senior bureaucrat turnover", "Including ministry switching in senior bureaucrat turnover",                          
                                "Including ministry switching \n and promotion/demotion in senior bureaucrat turnover")
modelsummary(OLS_alternative_DVs, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_alternative_DVs", "html"), 
             coef_rename = output_names,
             gof_map = gm,
             title = "OLS regression results: Different DV operationalizations",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

## Political buffer 
OLS_alternative_DVs_polbuff <- feols(c(dep.rad_turnover, dep.rad_turnover2, dep.rad_turnover3) ~ GOV_turnover * state_secretaries_NSDID +
                                       election_year + age_year + gender + temporary_employment + pre_unionsoppløsning | time+NSD_ID+decade, 
                                     cluster = "PersonID", df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
names(OLS_alternative_DVs_polbuff) <- c("Senior bureaucrat turnover", "Including ministry switching in senior bureaucrat turnover",                          
                                        "Including ministry switching \n and promotion/demotion in senior bureaucrat turnover")
modelsummary(OLS_alternative_DVs_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_alternative_DVs_polbuff", "html"), 
             coef_rename = output_names,
             gof_map = gm,
             title = "OLS regression results: Different DV operationalizations",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

## Calculating marginal effects and diff-in-diff
Coefplot_OLS_alternative_DVs <- model_marginal_effects(OLS_alternative_DVs) %>% mutate(model = DV_names)
Coefplot_OLS_alternative_DVs_polbuff <- model_marginal_effects(OLS_alternative_DVs_polbuff) %>% mutate(model = DV_names)

## Outputting Coefficient Plots
plotdata <- bind_rows(
  Coefplot_OLS_alternative_DVs, 
  Coefplot_OLS_alternative_DVs_polbuff
) %>% 
  pivot_longer(c(state_secretaries_NSDID, bureaucratic_buffer), values_drop_na = T) %>% 
  mutate(model = factor(model, levels = rev(c("Senior bureaucrat turnover",
                                              "Including ministry switching",
                                              "Including ministry switching \n and promotion/demotion"))),
         value = factor(value, levels = rev(c("Level-2",
                                              "Level-1",
                                              "Level-1 & No buffer",
                                              "Level-1 & Political buffer",
                                              "Difference-in-Differences"))))
print_plot(marginal_effect_plot(plotdata, models = T, facet = T), name = "coefplot_altDVs")


#### Figure 15 - Marginal effects Additional Moderator Variables ####
## Bureaucratic buffer 
OLS_robust_moderators <- feols(dep.rad_turnover ~ GOV_turnover * bureaucratic_buffer +
                                 election_year + age_year + gender + temporary_employment + pre_unionsoppløsning + 
                                 sw0(government_PMparty_share_posts_end_of_year, years_since_DEP_start, dep_end, multiple_L1, age_group_year2) | time+NSD_ID+decade, 
                               cluster = "PersonID", df)
names(OLS_robust_moderators) <- str_replace_all(c("Baseline", "government_PMparty_share_posts_end_of_year", "years_since_DEP_start", "dep_end", "multiple_L1", "age_group_year2"), output_names)
modelsummary(OLS_robust_moderators, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_robust_moderators", "html"), 
             coef_rename = output_names,
             gof_map = gm,
             title = "OLS regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

## Political buffer 
OLS_robust_moderators_polbuff <- feols(dep.rad_turnover ~ GOV_turnover * state_secretaries_NSDID +
                                         election_year + age_year + gender + temporary_employment + pre_unionsoppløsning + 
                                         sw0(government_PMparty_share_posts_end_of_year,years_since_DEP_start, dep_end, multiple_L1, age_group_year2) | time+NSD_ID+decade, 
                                       cluster = "PersonID", df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
names(OLS_robust_moderators_polbuff) <- str_replace_all(c("Baseline", "government_PMparty_share_posts_end_of_year", "years_since_DEP_start", "dep_end", "multiple_L1", "age_group_year2"), output_names)
modelsummary(OLS_robust_moderators_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}", 
             output = output("OLS_robust_moderators_polbuff", "html"), 
             coef_rename = output_names,
             gof_map = gm,
             title = "OLS regression results: DV = Turnover of senior bureaucrat",
             notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)

## Calculating marginal effects and diff-in-diff
Coefplot_OLS_robust_moderators <- model_marginal_effects(OLS_robust_moderators) %>% mutate(model = moderators_names)
Coefplot_OLS_robust_moderators_polbuff <- model_marginal_effects(OLS_robust_moderators_polbuff) %>% mutate(model = moderators_names)

## Outputting Coefficient Plots
Coefplot_moderators <- bind_rows(Coefplot_OLS_robust_moderators_polbuff,
                                 Coefplot_OLS_robust_moderators) %>% 
  pivot_longer(c(state_secretaries_NSDID, bureaucratic_buffer), values_drop_na = T) %>% 
  mutate(model = factor(model, levels = rev(str_replace_all(c("Baseline",
                                                              "government_PMparty_share_posts_end_of_year",
                                                              "years_since_DEP_start", "dep_end",
                                                              "multiple_L1", "age_group_year2"), output_names))))
print_plot(marginal_effect_plot(Coefplot_moderators, models = T), name = "Coefplot_moderators")
